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Abstract 



. We present results for the light quark masses obtained from a lattice QCD simulation 

' with Nf — 2 degenerate Wilson dynamical quark flavours. The sea quark masses of our 

lattice, of spacing a ~ 0.06 fm, are relatively heavy, i.e., they cover the range corresponding 
to 0.60 < Mp/Mv < 0.75. After implementing the non-perturbative RI-MOM method to 
rcnormalise quark masses, we obtain r7iJ^^(2GeV) = 4.3 ± 0.4tj!)'^ MeV, and m^^(2GeV) = 
lOliS^o^ MeV, which are about 15% larger than they would be if renormalised perturbatively. 
In addition, we show that the above results are compatible with those obtained in a quenched 
simulation with a similar lattice. 



1 Introduction 



An accurate determination of qTiark masses is highly important for both theoretical studies of 
flavour physics and particle physics phenomenology. Within the quenched approximation (Nj = 0), 
lattice QCD calculations reached a few percent accuracy, so that the error due to the use of 
quenched approximation became the main source of uncertainty. To examine the effects of the 
inclusion of the sea quarks several lattice results with either Nf = 2 or Nf = 3 dynamical quarks 
have recently been reported [1-9]. Many of them seem to indicate that the unquenching lowers the 
physical values of light quark masses. For example, the strange quark mass, which in the quenched 
approximation is about mf^(2GeV) ~ 100 MeV, becomes m^^^(2GeV) ~ 85 ^ 90 MeV when 
Nf = 2 dynamical quarks are included [2, 3], and even smaller with Nf = 3, namely m^^(2 GeV) ~ 
75^80 MeV [6,7]. 

The accuracy of the results obtained from simulations with Nf ^ 0, however, is still limited by 
several systematic uncertainties. In particular, precision quenched studies of light quark masses 
evidenced the importance of computing the mass renormalisation constants non-perturbatively. 
In almost all the studies with dynamical fcrmions performed so far, however, renormalisation has 
been made by means of one-loop boosted perturbation theory (BPT). Only very recently, non- 
perturbative renormalisation (NPR) has been implemented in a simulation with Nf = 2 by the 

QCDSF collaboration. The resulting strange quark mass value is mf^{2 GeV) = 119 it 9 MeV [8], 
thus larger by about 20% than the one previously reported in ref. [3], where the same lattice action 
and the same vector definition of the bare quark mass have been used, but with the renormalisation 
constant estimated perturbatively. Recently, also the Alpha collaboration reported their value for 
the NPR-ed strange quark mass, mY^{2 GeV) = 97 ± 22 MeV [9]. It appears that the results of 
both QCDSF and Alpha are consistent with their previous quenched estimates [10, 11]. 

In this paper we present our results for the light quark masses obtained on a 24^ X 48 lattice 
with Nf = 2 dynamical mass-degenerate quarks, in which we implement the RI-MOM NPR 
method of refs. [12, 13]. In the numerical simulation we used the Hybrid Monte Carlo algorithm 
(HMC) [14, 15], and chose to work with the Wilson plaquette gauge action and the Wilson quark 
action at /3 = 5.8. The corresponding lattice spacing is a ~ 0.06 fm, and the spatial extension of 
our lattice is L ~ 1.5 fm. Gauge configurations have been generated with four different values of 
the sea quark mass, for which the ratio of the pseudoscalar over vector meson masses lies in the 
range Mp/My — 0.60^ 0.75. Our final results for the average up/down and strange quark masses 
are 

mJS(2GeV) = 4.3(4)(+J-i) MeV / Nf = 2 

mp(2GeV) = 101(8)(1^^) MeV ^ RI-MOM 

where the first error is statistical and the second systematic, the latter dominated by discretisation 
lattice artifacts. Our central values refer to the quark mass definition based on the use of the axial 
Ward identity. The above results are larger than the ones we obtain by using one-loop BPT to 
compute the mass renormalisation, namely 

m^(2GeV) = 3.7(3)ej-^)MeV f Nf = 2 

mp(2 GeV) = 88(7) (1^0) MeV y BPT 

The latter values agree with the Nf =2 results reported in refs. [2, 3] in which BPT renormalisation 
has been used, while the values in eq. (1) agree with refs. [8, 9] in which NPR has been implemented. 
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The values given in eq. (1) can be also compared with our quenched estimate, 



mJf/(2GeV)=4.6(2)(tO-5)MeV / Nf = 

mf^(2GeV) = 106(2) (tj^ ) MeV \ RI-MOM 



(3) 



obtained on a lattice similar in size and resolution, and by employing the same Wilson quark 
action. We therefore conclude that, within our statistical and systematic accuracy, and with the 
sea quark masses used in our simulations, we do not observe any significant effect due to the 
presence of dynamical quarks. We plan to pursue this study by exploring lighter sea quarks. 

The remainder of this paper is organised as follows: in sec. 2 we discuss the details of our 
numerical simulations; in sec. 3 we present the results for the meson masses and the bare quark 
masses, directly accessed from the simulation; in sec. 4 wc provide the mass renormalisation factors, 
while in sec. 5 we determine the physical quark mass values; we briefly conclude in sec. 6. 

2 Simulation details 

In our numerical simulation we choose to work with the Wilson action 

Sw = Sg + Sq , (4) 
where Sg is the standard Wilson plaquette gauge action, 

^5 = ^E^^P = ^E^ [U.,,U,+ii,.ui_,,^^ul,] , (5) 
and Nf = 2 flavours of mass-degenerate dynamical quarks are included by using the Wilson action 



Sq — ^ ^ (IxDxyQy — ^ ^ Qx 
x,y x,y 



Sxy - l^'^ii^- 'llx)Ux,nSx+ii,y + (1 + 'ltl)Uli^dx,y+il\ 



Qy > (6) 



with K being the usual hopping parameter. The simulation parameters used in this work are 

summarised in table 1. To generate the ensembles of gauge field configurations we employ the 
HMC algorithm [14, 15] and implement the LL-SSOR preconditioning [16]. In solving the molecular 
dynamics equations we use the leapfrog integration scheme. Each trajectory of unit length (r = 1) 
is divided into iVgtep = 250 time-steps. The resulting St = r/iVstep = 4 x 10~^ appears to be 
small enough to avoid large energy differences along a trajectory, thus increasing the acceptance 
probability. Wc find an acceptance probability of about 85%. To check the reversibility along 
the trajectory wc verify that the quantity — 1|| remains unchanged at the relative level of 

10~^^. For the quark matrix inversion algorithm we used the BiConjugate Gradient Stabilized 
(BiCGStab) method [17]. The inversion rounding errors are of the order of 10~^. 

To select sufficiently dccorrclated configurations we studied the autocorrelation as a function 
of the number of trajectories t by which two configurations are separated. For a given observable 
A, the autocorrelation is evaluated as 

-, TMC—t 

C-it) = ^1- , (7) 



2 





/3 = 5.8 


/? = 5.6 




(J.i-.).");) U.iO-io (J.i-)4U U.i;)4i 


U.ioljU U.i-)/-.) l).i.)oU 


L '^ xT 


24^ X 48 


16^ X 48 


mimitcs/traj. 
Tmc 

^conf 
{P) 


10 12 16 21 
2250 2250 2250 2250 
50 50 50 50 
0.59268(3) 0.59312(2) 0.59336(3) 0.59341(3) 


5.3 13 28 
2250 2250 2250 
50 50 50 
0.56990(6) 0.57248(6) 0.57261(6) 


xT 


16^ X 48 




minutes/traj. 
Tmc 

^conf 
{P) 


3.3 4.4 4.4 5.6 
4500 4500 4500 4500 
100 100 100 100 
0.59283(3) 0.59314(3) 0.59340(3) 0.59345(3) 



Table 1: Run parameters for the simulations with Nf — 2. From the entire ensemble of generated HMC 
trajectories of unit length, Tmc, w^e select Nconf for the data analysis. The stopping condition in the 
quark matrix inversion (D^yGy — B^) is \\DG — < 10""'^^. (P) is the value of the average plaquette. 
"Minutes/traj. " refers to the time required to produce one trajectory on the APEmille machine (128 GFlop). 



where (^4) is the mean value of A computed on the total number of trajectories Tmc and Ag 
denotes the value measured along the trajectory s. After examining the aTitocorrclation of the 
plaquette and of the pseudoscalar two-point correlation functions, both shown in fig. 1, we made 
the conservative choice to select configurations separated by 45 trajectories. ^ The sea quark 
masses used in our simulations correspond to ratios of pseudoscalar over vector meson masses 
covering the range 0.60 < Mp/My < 0.75. In our data analysis we distinguish the sea (ks) from 
the valence (k,.) hopping parameters, where, for each value of the sea quark mass, Kj, can assume 
any of the values listed in table 1. The main results of this paper refer to the run made at /? = 5.8 
and L = 24 (i.e., L ~ 1.5 fm), whereas those obtained on the lattice with L = 16 and /? either 
5.8 (L 2± 1 fm) or 5.6 {L 2± 1.3 fm) are used to investigate finite volume and discretisation effects 
respectively. 

Finally, in order to assess the effects of quenching, we also made a quenched run (250 config- 
urations) at (3 = 6.2 and V = 24^ x 56 (L ~ 1.6 fm). To reproduce the values of the Nf = 2 
pseudoscalar meson masses, we chose in the quenched run G {0.1514,0.1517,0.1519,0.1524}. 

An important feature of the data analysis in the partially quenched case is the appropriate 
treatment of the statistical errors. While the results of computations of an observable O on 
configurations with different values of Kg are statistically independent, those computed on the 
gauge configurations with the same Kg but diflFerent Ky are not. To account for these correlation 
effects we use a combination of the bootstrap and jackknife statistical methods. At fixed Kg, our 
data are clustered into Nj^ subsets and the statistical errors on O are estimated by using the 
jackknife method. The resulting Nj^ values at a given Kg, 0^|, are then considered as a set of 
uniformly distributed numbers. In a next step one generates Ni, S> Njk bootstrap events by taking, 
for each Kg, one of the uniformly distributed O^^. The variance of O is then estimated by 

4 = {Njk - 1) X (^{0')n, - (O)ir,) , (8) 

^Tho autocorrelation of the pseudoscalar two-point function is made for the time separation between the sources 
fixed to 14. 
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PLAQUETTE AUTOCORRELATION PSEUDOSCALAR AUTOCORRELATION 




Figure 1: Autocorrelations of the plaquette C^{t) (left), and of the pseudoscalar correlation function C^^{t) 
(right), as a function of the number of trajectories t used to separate consecutive configurations. For C^^(t) 
the fluctuations around zero start at about t = 30, indicating that configurations separated by less than 30 
trajectories are correlated. The above plots refer to the run at (3 = 5.8 (V = 24"^ x 48, and Kg = 0.1538j. 

where {■■■)ni, stands for the average calculated over the Nf, bootstrap events, while the factor 
(Njk — 1) takes into account the jackknife correlation at fixed Hs- For a comparison, we also 
considered the procedure adopted in ref. [2], consisting in applying the jackknife method at one 
specific value of Kg, while keeping the data at the remaining three Kg fixed at their central values. 
As expected, after adding in quadrature the four variances calculated in this way, we find excellent 
agreement with the estimates obtained by using eq. (8). 

3 Masses of light mesons and quarks 

In this section we present the results for the light pseudoscalar and vector meson masses, as well as 
the bare quark masses, obtained from our main data-set (/3 = 5.8, L = 24). These quantities are 
extracted from the standard study of two-point correlation functions, Cjj{t), in which we consider 
both degenerate and non-degenerate valence quarks, with valence quarks either equal or different 
from the sea quarks. At large Euclidean time separations, t ^> 0, Cjj{t) is dominated by the 
lightest hadronic state and we fit it to the form: 

CjAt) = Y^{0\J{x,t)j\0m ^ ^(e-^^-^* + r?e-*^^(^-*)) , (9) 

X 

where rj is the time reversal {t ^ T — t) symmetry factor. ^ The pseudoscalar and vector me- 
son masses are extracted from Cpp{t) and Cvv{t) respectively, by fitting in t G [14,23]. The 
corresponding results are listed in table 2. These results are obtained by using local source opera- 
tors. We checked, however, that the Jacobi smearing procedure [18] leaves the masses unchanged, 
although the lowest hadronic state becomes isolated earlier. 

In the same table 2 we also give the results for the quark masses obtained by using the axial 

^77 — +1 for Cpp{t), Cvv{t), or CAA{t), whereas for CAp{t), t] = —1. In this notation P — qj^q, V = Vi = q'yiq, 
and A = Ao = q'yojsq. 
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Kg 






Mp 


My 




0.1535 


0.1535 


0.1535 


0.262(4) 


0.348(8) 


0.0333(6) 




0.1535 


0.1538 


0.252(5) 


0.341(8) 


0.0306(6) 




0.1535 


0.1540 


0.245(5) 


0.336(9) 


0.0289(6) 




0.1535 


0.1541 


0.241(5) 


0.334(9) 


0.0280(6) 




0.1538 


0.1538 


0.241(5) 


0.333(9) 


0.0280(6) 




0.1538 


0.1540 


0.233(5) 


0.329(9) 


0.0263(6) 




0.1538 


0.1541 


0.229(5) 


0.326(9) 


0.0254(6) 




0.1540 


0.1540 


0.226(5) 


0.324(9) 


0.0246(6) 




0.1540 


0.1541 


0.222(5) 


0.322(10) 


0.0237(6) 




0.1541 


0.1541 


0.218(5) 


0.319(10) 


0.0228(6) 


0.1538 


0.1535 


0.1535 


0.258(3) 


0.335(4) 


0.0314(2) 




0.1535 


0.1538 


0.247(3) 


0.328(4) 


0.0287(2) 




0.1535 


0.1540 


0.240(3) 


0.323(5) 


0.0269(2) 




0.1535 


0.1541 


0.236(4) 


0.321(5) 


0.0261(2) 




0.1538 


0.1538 


0.236(4) 


0.321(5) 


0.0260(2) 




0.1538 


0.1540 


0.228(4) 


0.316(6) 


0.0243(2) 




0.1538 


0.1541 


0.224(4) 


0.314(6) 


0.0234(2) 




0.1540 


0.1540 


0.220(4) 


0.311(6) 


0.0226(2) 




0.1540 


0.1541 


0.216(5) 


0.309(7) 


0.0217(2) 




0.1541 


0.1541 


0.212(5) 


0.306(7) 


0.0208(2) 


0.1540 


0.1535 


0.1535 


0.258(2) 


0.345(5) 


0.0303(3) 




0.1535 


0.1538 


0.247(2) 


0.338(6) 


0.0277(3) 




0.1535 


0.1540 


0.240(3) 


0.333(7) 


0.0259(3) 




0.1535 


0.1541 


0.236(3) 


0.331(7) 


0.0250(3) 




0.1538 


0.1538 


0.236(3) 


0.331(7) 


0.0250(3) 




0.1538 


0.1540 


0.228(3) 


0.326(8) 


0.0233(3) 




0.1538 


0.1541 


0.225(3) 


0.324(8) 


0.0224(3) 




0.1540 


0.1540 


0.221(3) 


0.321(9) 


0.0215(4) 




0.1540 


0.1541 


0.217(3) 


0.319(9) 


0.0207(4) 




0.1541 


0.1541 


0.212(3) 


0.316(10) 


0.0198(4) 


0.1541 


0.1535 


0.1535 


0.234(6) 


0.318(10) 


0.0286(3) 




0.1535 


0.1538 


0.222(6) 


0.311(11) 


0.0259(3) 




0.1535 


0.1540 


0.213(6) 


0.307(11) 


0.0242(3) 




0.1535 


0.1541 


0.209(6) 


0.304(11) 


0.0233(3) 




0.1538 


0.1538 


0.209(6) 


0.304(11) 


0.0233(3) 




0.1538 


0.1540 


0.200(7) 


0.300(12) 


0.0215(3) 




0.1538 


0.1541 


0.196(7) 


0.298(12) 


0.0206(3) 




0.1540 


0.1540 


0.191(7) 


0.295(12) 


0.0198(3) 




0.1540 


0.1541 


0.187(7) 


0.293(13) 


0.0189(3) 




0.1541 


0.1541 


0.182(7) 


0.291(13) 


0.0180(4) 



Table 2: Pseudoscalar and vector meson masses, and the bare quark masses m^i2^{a) = ^[m^i(a) + 

771^,2(0)]^^^ obtained by using the axial Ward identity method (cf. cq. (10)). Rcsjdts, in lattice units, refer 
to the simulation with f3 = 5.8 and V = 2^ x 48. Those obtained at P = 5.6 are tabulated in appendix 1. 
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Figure 2: Data-points that illustrate the functional dependencies discussed in eqs. (12) and (16) 



Ward identity (AWI) definition. More specifically, with Aq{x) = qvi7o75Qv2, and P = qvi^5qv2, 



f»0 



(j;P(x)pt(o)) 



[m„i(a) + m„2(a)] 



AWI 



(10) 



where we fit in the same interval as before, i.e., t S [14,23]. Alternatively, one can use the vector 
Ward identity (VWI) and relate the bare quark mass to the Wilson hopping parameter as: 



m. 



VWI 



(a) 



1 



1 



(11) 



with K„ being either (cf. table 2), and Kcr{K.s) is the critical hopping parameter for 

a given sea quark mass. Notice that both the AWI and VWI definitions of the quark mass do 
depend on the values of the sea quark masses. In the axial case, this dependence is provided by the 
ratio of two-point correlation functions of eq. (10). In the VWI definition, the dependence on the 
sea quark mass is given by the value of the critical hopping parameter Kcr(^s) computed at fixed 
Kg. We also notice that the sea quark dependence of the critical hopping parameter is numerically 
equivalent to the observation made in ref. [8] where such a dependence is expressed in terms of the 
renormalisation constants of the singlet and non-singlet scalar density. 



3.1 Lattice spacing and Kcri^g) 

In order to fix the value of the lattice spacing and Hct{i^s), we inspect the dependence of the 
pseudoscalar and vector meson masses on the valence and sea quark masses. The left panel of 
fig. 2 shows the dependence of the vector meson masses on the squared pseudoscalar ones, which 
suggests a simple linear fit to the form 

Mv{Ks, I^Vi,I^V2) = Pl+ P2X Mp{Ks,Ky^,K.V2) + P3X Mp{K.s, Kg, Kg) , (12) 

in an obvious notation, and with capital letters used to denote that the meson masses are given 
in lattice units. Such a fit yields P3 = 0.1(3), i.e., the sea quark dependence is very weak. We 
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can therefore neglect it and set P3 = 0. In this way we obtain Pi = 0.25(1) and P2 = 1.4(2). We 
then determine the lattice spacing from this fit and by using the method of the "physical lattice 
planes" [19]: at the intersection of our data with the physical value for the ratio rriK/mK* = 0.554, 
we impose that Mk* obtained on the lattice is Mk* = amF^»^ , thus obtaining 

a-^^, =3.2(1) GeV. (13) 

We also estimate the value of the lattice spacing by studying the static quark potential. At each 
value of the sea quark mass we extract the following values of the Sommer parameter Rq = ro/a [20] : 

(i?o)«, = {7.52(17)0.1535, 7.70(9)0.1538, 7.78(19)0.1540, 8.10(20)o.i54i}. (14) 

Details on this extraction can be found in appendix 2. Assuming the observed dependence of Rq 
on the sea quark mass to be physical, as the lattice spacing a is supposed to be independent on 
the sea quark mass, we linearly extrapolate the above results to the chiral limit [i?o = a + P x 
Mp{Ks, Kg, Kg)], ending up with Rq = 8.6(4). After setting ro = 0.5 fm, we find 

a-' = 3.4(2) GeV , (15) 

in good agreement with the value given in eq. (13). 

Finally, for m^^^ in eq. (11) we also need Kcr{Ks)- Inspired by partially quenched ChPT [21], 
but neglecting the chiral logs since our dynamical quarks are not light enough, we fit our data to 

mUi^s, k^„k^,) = Qi mX^'[l + Q2 mX^' + Q3 m^^'] , (16) 

and obtain Q2 = —0.2(25) and Q3 = 0.4(14). In other words our data are not sensitive to the 
quadratic corrections proportional to Q2 and Q3, as it can be seen from the right plot in fig. 2. 
Very similar feature is observed if instead of 111^^^ we use m^^^ , and therefore in the following 
we set Q2 = Q3 = 0. Prom the resulting fit to eq. (16), we get Qi = 1.70(3) and 

(«cr)«. = {0.15544(7)0.1535, 0.15537(6)0.1538, 0.15537(5)0.1540, 0.15503(8)o.i54i} • (17) 



4 Quark mass renormalisation 

Before discussing the strategy to identify the physical strange and the average up/down quark 
masses from the results obtained by using eqs. (10,11), we need to determine the corresponding mul- 
tiplicative mass renormalisation constants, Z:^^{aiJ,) = ZA/Zp{an), and Z^^^{aiJ,) = 1/Zs{an). 
For completeness, we also present in this section the results for Zy, Zt and the quark field RC Zq. 

We use the RI-MOM method [12] which we explained in great detail in a previous publica- 
tion [13]. We adopt the same notation as in ref. [13] and compute the renormalisation group 
invariant combinations 

Zoiaiio) = Coiafio) Z§^^ = Co{aiio) {Zo{aii)/Co{aii)) , (18) 

for the scale dependent bilinear operators, O = S,P,T. The evolution functions Co{a^i), which 
arc known in the RI-MOM scheme at the N^LO for Zs and Zp [22] and at the N^LO for Zt [23], 
explicitly cancel the scale dependence of the RCs, Zo{aiJ,), at the considered order in continuum 
perturbation theory. 

In fig. 3 wc plot Zy^A, SlS well as Zs^p^riO'fJ'o), where in eq. (18) wc choose afiQ = 1- All of them 
are supposed to be flat in {a/j,)^, where /x is the initial renormalisation scale introduced in eq. (18), 
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0.85 




0.75 



0.65 



0.55 



0.0 



0.5 



1.0 



1.5 



2.0 



Figure 3: Renormalisation constants in the RI-MOM scheme at the scale afiQ = 1, computed according to 
eq. (18), as a function of the initial renormalisation scale (afi)'^ . Discretisation effects oc (a/i)^ are eliminated 
by extrapolating the data from 1 < (a/i)^ < 2, to (a/i)^ = 0. Illustration is provided for the data with 



modulo discretisation effects. In fact, however, we see that these effects are pronounced in the 
case of Zsia^o), which we correct for by hnearly extrapolating from 1 < (a//)^ < 2, to (a//)^ = 0. 
Furthermore, for the RI-MOM scheme to be mass independent we extrapolate the mild sea quark 
dependence hnearly to the chiral limit. The results are reported in table 3. By confronting the 
results for the RCs obtained at /3 = 5.8 and L = 24 with those obtained at the same /? but with 
L = 16, we find that they differ by less than 5%. Given the smallness of the physical lattice size at 
L = 16 (~ 1 fm), we expect finite volume effects on our estimate of the RCs to be negligible with 
respect to discretization errors, which are evaluated to be of the order of 10% (see subsec. 5.1). 
This is expected, since the RCs encode the short distance physics while the volume finiteness affects 
the long distance physics. For an easier comparison of the non-perturbatively and perturbatively 
estimated RCs, in table 3 we also give the RCs evaluated by means of 1-loop BPT [24]. These are 
obtained by replacing in the one-loop expressions of ref. [25] the bare lattice coupling by a boosted 
one, 5o ~^ 5^- distinguish the so-called "naive" boosting 



with (P) being the average plaquette, and the so called "tadpole improved" coupling, which is 
defined by inverting the perturbative series of the logarithm of the plaquette, namely ^ 



The corresponding results for the RCs are then referred to as N-BPT and TI-BPT respectively. 
From the results in table 3 we see that the perturbative estimates of Zs^p are larger than the non- 
perturbative ones. In other words the non-perturbatively computed .^^^^ = Za/Zp (Z^^^ = 
1/Zs), at (3 = 5.8, is larger by about 15% (6%) than its BPT counterpart. These differences 
directly propagate to the lattice determination of the renormalised quark masses. As we shall see 

^For a more extensive discussion about various forms of boosting, see ref. [26]. 



K, = 0.1538. 



f = 91/ (P) 



(19) 




g'(3.40/a) 
(4vr) 



(20) 



8 





Kg 




Zv 


Za 


Zs 


Zp 




5.8 (Nf = 2) 


0.1535 


0.81(1) 


0.70(1) 


0.80(2) 


0.74(1") 


0.54(1) 


0.78(2) 




0.1538 


0.80(1) 


0.70(1) 


0.79(2) 


0.73(2) 


0.55(1) 


0.77(2) 




0.1540 


0.78(1) 


0.67(1) 


0.77(2) 


0.72(1) 


0.52(1) 


0.75(2) 




0.1541 


0.79(1) 


0.69(1) 


0.78(2) 


0.71(2) 


0.52(1) 


0.77(2) 




^cr ( ^cr ) 


0.76(1) 


0.66(2) 


0.76(1) 


0.67(3) 


0.50(2) 


0.74(4) 




N-BPT 


0.75 


0.70 


0.77 


0.75 


0.61 


0.75 




TI-BPT 


0.72 


0.65 


0.73 


0.71 


0.55 


0.71 


5.6 {Nf = 2) 


0.1560 


0.78(2) 


0.63(4) 


0.78(3) 


0.67(2) 


0.46(1) 


0.76(5) 




0.1575 


0.78(2) 


0.64(4) 


0.77(2) 


0.65(2) 


0.46(1) 


0.77(3) 




0.1580 


0.78(1) 


0.66(1) 


0.78(4) 


0.69(2) 


0.50(1) 


0.78(3) 




l^cr{l^cr) 


0.78(1) 


0.66(1) 


0.78(1) 


0.67(1) 


0.47(1) 


0.78(2) 




N-BPT 


0.74 


0.67 


0.75 


0.73 


0.58 


0.73 




TI-BPT 


0.69 


0.62 


0.71 


0.69 


0.51 


0.68 


6.2 {Nf = 0) 




0.82(1) 


0.71(1) 


0.80(1) 


0.71(1) 


0.54(1) 


0.79(3) 




N-BPT 


0.78 


0.73 


0.79 


0.77 


0.65 


0.77 




TI-BPT 


0.72 


0.64 


0.73 


0.71 


0.55 


0.71 



Table 3: Values of the renormalisation constants in the RI-MOM scheme at the scale ^ ~ 1/a for the 
two partially quenched simulations at (3 = 5.8 (a~^ = 3.2(1) GeVj and 5.6 (a~^ = 2.4(2) GcVj and for 
the quenched simulation at (3 ~ 6.2 (a^^ = 3.0(1) GcVj. The quoted errors arc statistical only. The 
non-perturbatively determined RCs are confronted to their perturbative value (BPT), where "TF and "N" 
stand for the two types of boosting the bare lattice coupling as discussed in the text. Definition of Zq is the 
same one as in eq. (5) of ref. [13]. 



in the following this effect explains the differences between our non-perturbatively renormalised 
results and those obtained by the CP-PACS and JLQCD Collaborations [2,3], where the TI-BPT 
has been used. Furthermore, our results agree with those obtained in refs. [8,9], in which the 
non-perturbative renormalisation has been implemented. 

5 Physical light quark masses 

In order to determine the physical values of the light quark masses, we need to investigate the 
dependence of the pseudoscalar meson masses squared on the valence and sea quark masses. Prom 
the left plot in fig. 4 we see that: (i) the dependence on the valence quark mass is linear, (ii) the 
slopes do not depend on Ks, and (iii) the intercepts do depend on Kg and they are non-zero even 
when the sea quark mass is sent to zero. We can therefore fit to the form 

MUi^g, K,„K,,) = Co{ks) + Ci m^'^'il + C2 mj"^' + C3 m^"^'] , (21) 

where the valence and sea AWI quark masses are defined in eq. (10). In this notation m^^^ = 
m^^\Ks, K^^, Ki,,^), and m^^^ = rn"^^^(K5, k^,, k^). Obviously, the fit forms (16) and (21) differ 
not only in the use of differently defined quark masses (AWI vs. VWI), but also in the presence 
of the constant term Cq{ks) in eq. (21). This term accounts for the 0(a)-discretisation effects 
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Figure 4: Dependence of the pseudoscalar meson masses on the valence quark masses m^^^ (left) and of 
the ratios [Mp — Co)/m^^-' (cf. eq. (21)) on the sea quark masses mf^^ . The results refer to the case of 
the partially quenched simulations at (3 = 5.8. 



present in either the pseudoscalar meson masses or the AWI quark masses or in both. Clearly, 
discretisation effects in the pseudoscalar meson masses also affect the fit to eq. (16), i.e., the 
determination of the critical hopping parameters Kcr{i^s)- 

As before, we find that the values of the coefficients of the quadratic terms, C2 and C3, in 
eq. (21) are consistent with zero [C2 = 0.6(18), and C3 = 0.08(63)], which can be also seen from 
the right plot in fig. 4. Therefore we may set C2 = C3 = 0, and from the fit to eq. (21) we obtain 
Ci = 2.05(3), and 

Co(k,) = {0.0007(18)0.1535, 0.0023(20)0.1538, 0.0046(15)o.i540, -0.0039(25)o.i54i} • (22) 

Our partially quenched estimates of the average up/down {rriud) and of the strange {nis) quark 
masses are made by substituting on the l.h.s. of eqs. (16) and (21) the physical pion and kaon 
masses and using the lattice spacing value given in either eq. (13) or eq. (15). In other words, to 
get the quark masses in physical units, we solve 

{miy^^"' =QiX X ruud , 

^^2^Yhys. ^ ^ ^ ^-1 ^ ^ ^ (23) 

in the VWI case, and the analogous expressions, with the coefficient Qi substituted by Ci, in 
the AWI case. Finally the bare quark masses are multiplied by the corresponding RI-MOM mass 
renormalisation constants computed at ^0 = l/o, already listed in table 3. 

The conversion to the reference MS scheme, at // = 2 GeV, is made by using the perturbative 
expressions known up to N^LO accuracy [22,27], and A^^"^^ = 0.245(16)(20) MeV [28]. We 
finally obtain 

VWI : m^(2 GeV) = 5.1(4) MeV , mf^{2 GeV) = 120(7) MeV , 

AWI : m^(2 GeV) = 4.3(4) MeV , mp(2 GeV) = 101(8) MeV . (24) 

Notice that our VWI result is in very good agreement with the value m^^(2 GeV) = 119(5)(8) MeV 
obtained by the QCDSF-UKQCD collaboration who uses the VWI quark mass definition and the 
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RI-MOM non-perturbative renormalisation [8]. Notice also that our AWI result agrees with the 
value obtained by the Alpha collaboration, mf^(2GeV) = 97(22) MeV [9]. 

5.1 Systematic uncertainties 

o Discretisation I : Prom quenched studies we learned that a discrepancy between the quark 
mass values obtained by using two a priori equivalent definitions, namely the VWI and 
AWI definitions, is due to discretisation errors. Such a difference has also been observed in 
other lattice studies with Nf = 2 [2,3]. In particular, in ref. [2] this difference is shown to 
diminish as the lattice spacing is reduced so that, in the continuum limit, the quark masses 
defined in the two ways converge to the same (unique) value. From that study it became 
clear that discretisation errors are significantly larger in the VWI determination of the quark 
masses which is why the JLQCD collaboration choose to quote their AWI masses as final 
results, while the difference between the VWI and AWI results is included in the asymmetric 
systematic error. We follow the same procedure here and we obtain 

m^(2 GeV) = 4.3(4)(+°-8) MeV , mp(2 GeV) = 101(8)(1^9) MeV . (25) 

o Discretisation II : The results reported in eq. (24) refer to the lattice with coupling (3 = 5.8 
(and volume 24^ x 48). Prom our data produced at /? = 5.6 (see table 1) we obtain 

m^(2 GeV) = 4.2(3)(+°-^) MeV , mf^(2 GeV) = 101(6)(1J^) MeV , (26) 

thus indistinguishable from those obtained at (3 = 5.8. This makes us more confident that 
the discretisation error attributed to our results in eq. (25) is realistic. 

o Lattice spacing: To obtain the results quoted in eq. (25) we used the lattice spacing (13) 

fixed by mF^T' . These results would become larger by 6% if we used (15). We will add 
this difference to our final error budget. 

o Renormalisation constants : As discussed in the text the determination of the mass renor- 
malisation factors is sensitive to discretisation errors. They are especially important when 
the unimproved Wilson quark action is used, which is the case with our study. To check for 
these effects we determined Zp/Zs by using the hadronic Ward identities (hWI) [29] (see 
eqs. (16,17) of ref. [13]), obtaining 

(Zp/Z5)^^' = 0.67(1). (27) 

Compared to the results in table 3, where aX (3 = 5.8 we find {Z p / Z s f-^'^^^ = 0.75(1), this 
result is about 11% smaller. A similar effect is seen at /3 = 5.6. Wc observe that, either one 
attributes this effect to Zp or to Z5, {Zp/Zs)^^^^ always goes in the direction of decreasing 
the AWI- VWI mass difference. Therefore, this 11% uncertainty, being already enclosed in the 
19% systematic error due to the AWI- VWI mass difference, will not be added. We finally 
notice that the result obtained for Zp/Zg by using the numerical stochastic perturbation 
theory to 4-loops, i.e., {Zp/Zs)^^^'^ = 0.77(1) [30], is more consistent with our larger NPR- 
ed (i.e., RI-MOM) value. 

o Finite volume : Our data at /? = 5.8 and 16^ x 48 (see appendix 1) indicate large finite volume 
effects, when compared to our main data-set with /3 = 5.8 and 24^ x 48. As an illustration 
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Figure 5: Pseudoscalar and vector effective mass plots at two different volumes at P = 5.8 and with 
Ks = 0.1538. Valence quark masses are equal to that of the sea quark. The physical volumes correspond to 
(1.5 fm)'^ and (1 fm)'^ respectively. 



we show in fig. 5 the effective mass plots obtained in the two volumes at the same value of 
the sea quark mass. Such a sensitivity to the finiteness of the lattice box has already been 
observed in ref. [31], where it was shown that, for the range of sea quark masses used in our 
study, finite volume effects are large at L ~ 1 fm but become insignificant when working on 
lattices with L > 1.5 fm. We will rely on that conclusion and will assume that finite volume 
effects are negligible with respect to our statistical and other sources of systematic errors. 
For a more refined estimate of these effects calculations on lattices with larger volumes are 
required. 

We add the above systematic errors algebraically and, as a final result, we quote 

m^(2 GeV) = 4.3 ± 0.4+^'^ MeV f Nf = 2, /3 = 5.8 

mp(2 GeV) = 101 ± 8+o^ MeV y RI-MOM 

which are the values already given in the abstract and the introduction of the present paper. 




5.2 Impact of non-perturbative renormalisation and quenching effects 

Before concluding this section we should compare our main results for the light quark masses 
quoted in eq. (28) with those obtained by evaluating the quark mass RCs using one-loop tadpole 
improved BPT: 

m^(2GeV) = 3.7±0.3lo^ MeV / iV^ = 2 , /? = 5.8 

mp(2 GeV) = 88 ± 7tf MeV y TI-BPT 

These results are in good agreement with those available in the literature obtained from lattice stu- 
dies with Nf = 2 in which renormalisation is implemented perturbatively, namely m^^(2GeV) = 
88(1^) MeV [2], and mf^(2GeV) = 84.5(^1^) MeV [3]. Therefore, with Wilson-like dynamical 
fermions, the non-perturbative renormalisation increases the values of the quark masses. 
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Finally, when compared with our results obtained in the quenched simulation at /3 = 6.2, 



on a lattice similar in size and resolution and by using the same (unimproved) Wilson quark action, 
we see no evidence for any effect that could be attributed to the presence of the dynamical quarks. 
For a clearer conclusion concerning this point one should work with lighter sea quark masses which 
is one of the main goals of our future simulations. 

6 Conclusions 

In this paper we have presented our results for the light quark masses obtained from a numerical 
simulation of QCD on the lattice with Nf = 2 degenerate dynamical quarks, with masses covering 
the range 3/4 < mq/nis^^^ ^3/2 (with respect to the physical strange quark mass). Our main 
values for the strange and for the average up/down quark masses are obtained on a lattice with 
spacing a ~ 0.06 fm and spatial volume (1.5 fm)^, and by using the axial Ward identity. An 
important feature of our calculation is that we renormalise the quark masses non-perturbatively 
in the RI-MOM scheme, which is then converted to the MS scheme by using the known 4-loop 
perturbative formulae. We show that 

— within our statistical accuracy and with relatively heavy dynamical quarks, our unquenched 
results are fully consistent with the quenched ones; ^ 

— the non-perturbative renormalisation leads to resulting quark masses larger than those renor- 
malised by using 1-loop perturbation theory (even if tadpole improved); 

— the systematic errors are dominated by lattice discretisation artifacts, which can be cured by 
implementing C'(a)-improvement and/or working with several small lattice spacings followed 
by an extrapolation to the continuum limit. 

The finiteness of the lattice volume is expected not to be an important source of systematic errors 
for the sea quark masses explored in our simulation. Further decrease of the dynamical quark 
masses would necessitate, however, to work with larger physical lattice volumes. Our final results 
are given in eq. (28). 
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mJf/(2 GeV) = 4.6 ± 0.2l°-^ MeV 

mp(2 GeV) = 106 ± 2+^^ MeV 




(30) 
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Appendix 1: Meson masses from two simulations with larger cou- 
pling {/3 = 5.6) and smaller volume (16^ x 48) 



In this appendix wc provide the reader with two sets of results. Wc first give our values equivalent 
to the ones presented in table 2 but for (3 = 5.6. Then, for the reader to better appreciate the size 
of finite volume effects at /3 = 5.8, we also tabulate the pseudoscalar and vector meson masses on 
the spatial volumes 16^ and 24? in the fully unquenched situations, i.e., with 



Kg 






Mp 


Mv 




0.1560 


0.1560 


0.1560 


0.441(5) 


0.541(6) 


0.0669(4) 




0.1560 


0.1575 


0.398(6) 


0.513(8) 


0.0548(5) 




0.1560 


0.1580 


0.383(7) 


0.503(8) 


0.0509(5) 




0.1575 


0.1575 


0.351(7) 


0.482(9) 


0.0430(5) 




0.1575 


0.1580 


0.334(8) 


0.471(10) 


0.0392(6) 




0.1580 


0.1580 


0.317(9) 


0.459(10) 


0.0355(6) 


0.1575 


0.1560 


0.1560 


0.379(3) 


0.456(8) 


0.0515(4) 




0.1560 


0.1575 


0.333(4) 


0.417(13) 


0.0394(4) 




0.1560 


0.1580 


0.316(4) 


0.403(16) 


0.0355(4) 




0.1575 


0.1575 


0.280(4) 


0.374(21) 


0.0276(3) 




0.1575 


0.1580 


0.261(4) 


0.357(27) 


0.0238(3) 




0.1580 


0.1580 


0.241(4) 


0.340(35) 


0.0200(3) 


0.1580 


0.1560 


0.1560 


0.368(3) 


0.468(11) 


0.0466(4) 




0.1560 


0.1575 


0.324(5) 


0.443(15) 


0.0348(4) 




0.1560 


0.1580 


0.308(6) 


0.436(18) 


0.0311(4) 




0.1575 


0.1575 


0.271(7) 


0.419(21) 


0.0234(5) 




0.1575 


0.1580 


0.250(9) 


0.419(27) 


0.0199(5) 




0.1580 


0.1580 


0.227(11) 


0.424(37) 


0.0164(5) 



Table 4: Values of pseudoscalar and vector meson masses, as well as of the AWI quark masses m^^^{a) = 
\ \mvi{a) + m„2(a)]^^^ obtained by using cq. (10), as obtained from the simulation at /3 = 5.6. All results 
are given in lattice units. 



Kg Ky 


M^ M^^'^ Ap 


Mr M^''^ Ay 


0.1535 
0.1538 
0.1540 
0.1541 


0.262(4) 0.299(6) 0.14(4) 
0.236(4) 0.272(10) 0.15(4) 
0.221(3) 0.257(8) 0.17(4) 
0.182(7) 0.251(9) 0.38(4) 


0.348(8) 0.402(13) 0.16(5) 
0.321(5) 0.356(17) 0.11(6) 
0.321(9) 0.347(23) 0.08(7) 
0.291(13) 0.328(31) 0.13(8) 



Table 5: Pseudoscalar and vector meson masses, in lattice units, calculated at (3 = 5.8 on the lattices 
24^ X 48 and 16^ x 48, respectively. Finite volume effects are quantified by Apy = {M^p^^ - M^p^) / M'^pl^ . 
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Appendix 2: The scale parameter Rq = ro/a 



To determine the Sommer scale parameter Rq [20] we computed the static quark potential V{R) 
which is extracted from the time dependence of the rectangular Wilson loops W{R,t), where R is 
the spatial separation between the static charges, i.e., 



WiRJ) = C(fl)e-"<«)' ^ V(R) = l^^VMRJ) ^ log {^^^) 



(31) 



Large time separations are reached after using the so-called HYP (hypercubic blocking) procedure 
proposed in ref. [33] which helps in improving the signal for the Wilson loops. In this way we find 
a stable signal for V{R) for tgt € [8,13]. At intermediate spatial separations between the static 
charges, we then fit the data (shown in fig. 6) to the form 



V{R) = Vo + aR--- g6V{R) , 



(32) 



where we also account for SV{R), the perturbatively estimated discretisation correction to the 
Coulomb potential [34]. The Sommer scale Rq is defined as 



R 



2dViR)\ 



dR 



) 



1.65. 



(33) 



R=Ro 



The fitting window consistent with the form (32) is found for R G [2, 7] after which the statistical 
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Figure 6: Static quark potential obtained from the simulation with Nf = 2 at (3 = 5.8 on the lattice 
X 48. 



quality of the data for V{R) rapidly deteriorates. From the fit of our data to the form (32) in that 



15 



window, we obtain 



(3 = 5.8 : Ks = {0.1535, 0.1538, 0.1540, 0.1541} , 

Ro = {7.52(17), 7.70(9), 7.78(19), 8.10(20)}, 

= {0.155(4), 0.151(2), 0.150(4), 0.144(4)}, (34) 

e = {0.289(10), 0.302(5), 0.294(7), 0.288(13)}, 

Vo = {0.260(7), 0.268(3), 0.264(5), 0.262(9)}, 

g = {0.033(14), 0.036(9), 0.035(7), 0.015(15)}, 

ordered as the values of Kg. Our results for Rq obtained at the same value of (3 but on the smaller 
volume (16^ x 48) agree with those written in eq. (34), within the statistical errors. 

At /? = 5.6, we find stability for tf^ G [3, 5], and then from the fit to eq. (32) in G [2, 7], we 
have 



/3 = 5.6: Ks = {0.1560,0.1575,0.1580}, 

Ro = {5.15(5), 5.72(10), 5.98(8)}, 

Va = {0.225(3), 0.202(4), 0.195(3)}, (35) 

e = {0.305(11), 0.317(8), 0.308(6)}, 

Vo = {0.262(9), 0.279(7), 0.275(5)}, 

g = {0.009(11), 0.022(7), 0.020(4)}, 

in good agreement with ref. [31]. 
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